Thermal dissociation of dipositronium: path integral Monte Carlo approach 
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Path integral Monte Carlo simulation of the dipositronium "molecule" Ps2 reveals its surprising 
thermal instability. Although, the binding energy is ~ 0.4 eV, due to the strong temperature 
dependence of its free energy Ps2 dissociates, or does not form, above ~ 1000 K, except for high 
densities where a small fraction of molecules are in equilibrium with Ps atoms. This prediction is 
consistent with the recently reported first observation of stable Ps2 molecules by Cassidy & Mills 
Jr., Nature 449, 195 (07), and Phys.Rev.Lett. 100, 013401 (08); at temperatures below 1000 K. 
The relatively sharp transition from molecular to atomic equilibrium, that we find, remains to be 
experimentally verified. To shed light on the origin of the large entropy factor in free energy we 
analyze the nature of interatomic interactions of these strongly correlated quantum particles. The 
conventional diatomic potential curve is given by the van der Waals interaction at large distances, 
but due to the correlations and high derealization of constituent particles the concept of potential 
curve becomes ambiguous at short atomic distances. 

PACS numbers: 31.15.xk, 36.10.Dr, 31.15.ae 



Dipositronium or positronium molecule, Ps2, is a 
four-body system consisting of two electrons and two 
positrons. The dynamical stability of dipositronium was 
established in 1947 by Hylleraas and Ore [lj]. How- 
ever, the molecule was not observed experimentally un- 
til recently Q, even though a lot of knowlegde had 
been provided by a number of theoretical studies, see 
Refs. [H, 0, H, E3, 0, @, 0] and references therein. In addi- 
tion to the fundamental issues of physics, Ps2 is of interest 
also in astrophysical applications and solid state physics 

EES. 

In laboratory conditions, Ps2 formation has recently 
been observed resulting from implantation of intense 
pulses of positrons into porous silica films 0, EH • 

The positronium molecule, with all the four parti- 
cles of the same mass, sets challenges to modeling, 
since quantum calculations are to be performed fully 
non-adiabatically This, however, can be realized 

with quantum Monte Carlo (QMC) methods [3, 
EH- It should be pointed out that also for other sys- 
tems, approaches that are not restricted by the Born- 
Oppenhcimer or other adiabatic approximations are gain- 
ing more attention 0, El El EE E3, EH, [H, El • 

Among the QMC methods the path integral Monte 
Carlo (PIMC) offers a finite temperature approach to- 
gether with a transparent tool to trace the correlations 
between the particles involved. Though computationally 
challenging, with the carefully chosen approximations 
PIMC is capable of treating low-dimensional systems, 
such as small molecules or clusters accurately enough for 
good quantum statistics for a finite temperature mixed 

state g, m, m E3, m si. 

In this study, using PIMC we evaluate the density ma- 
trix of the full four-body quantum statistics in tempera- 
ture dependent stationary states. Thus, the temperature 
dependent distributions of structures and energetics of 
Ps2 are established. The main focus here is to find the 



preferred configuration of the four-body system at each 
temperature — Ps2 molecule or two Ps atoms. 

According to the Feynman formulation of the statis- 
tical quantum mechanics [30| the partition function for 
interacting distinguishable particles is given by the trace 
of the density matrix p((3) 
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where S is the action, (3 — 1/ksT, r = (3/M and M is 
called the Trotter number (Rm = Ro)- In present simu- 
lations, we use the pair approximation of the action and 
matrix squaring for evaluation of the Coulomb interac- 
tions [25l.l3ll]. Sampling of the paths in the configuration 
space is carried out using the Metropolis algorithm [32| 
with the bisection moves 33] . The Coulomb potential en- 
ergy is obtained as an expectation value from sampling 
and the kinetic energy is calculated using the virial esti- 
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mator 

The error estimate for the PIMC scheme is comm only 
given in powers of the imaginary time time-step r [25(. 
Therefore, in order to determine comparable thermal ef- 
fects on the system we have carried out the simulations 
with similar sized time-steps regardless of the tempera- 
ture. This way the temperature dependent properties can 
be compared avoiding temperature dependent systematic 
errors. The standard error of the mean (SEM) with two- 
sigma limits is used to indicate the statistical uncertainty, 
where relevant. The average of the chosen time-step is 
(t) « 0.0146-Ejj 1 , where En denotes the atomic unit of 
energy, Hartree (~ 27.2 eV). The other atomic unit we 
use here is Bohr radius for the length, a (~ 0.529 A). 

The total energy of positronium "atom" Ps is —0.25 
at K and the binding energy of the molecule PS2 is 
0.0160 (« 0.435 eV) [|. We find these values as zero 
Kelvin extrapolates from our simulations at low temper- 



2 



atures. We point out that with PIMC we evaluate en- 
ergetics as statistical expectation values from sampling 
with less accuracy than that from conventional solutions 
of wavefunctions and the zero Kelvin data we obtain as 
extrapolates, only. 

In Fig. [1] we present the "apparent dissociation energy" 
of PS2 at several different temperatures. In each tempera- 
ture this is the negative total energy of the molecule with 
respect to two atoms as D T = -[(E^ t 2 ) T -2(E^ t ) T ]- At 
T < 900 K we find for the average over shown temper- 
atures Dt — 0. Of 54 (5), which is very close to the dis- 
sociation energy at zero Kelvin, Dq. However, at higher 
temperatures the apparent dissociation energy vanishes, 
because (E^)t and 2(_E t ^ s t )T become the same. This is 
because of molecular dissociation, or to be more exact, 
the two atoms do not bind in our equilibrium state sim- 
ulation at T > 900 K and the predominant configuration 
is that of two separate positronium atoms. 

Simulations in a well-defined Ps density are time con- 
suming, and therefore, these kind of studies have been 
carried out at the transition region around f 000 K, only. 
Using the periodic boundary conditions and the cubic 
supercells of sizes from (300ao) 3 to (50ao) 3 with two 
Ps atoms we have simulated three densities from 0.5 to 
100 x 10 24 m _1 , respectively. We see that with increasing 
density the equilibrium shifts to the molecular direction 
making the transition smoother and raising it to higher 
temperatures compared to the more sharp low density 
limit. 

For completeness we should point out that in equi- 
librium at any finite temperature the zero density limit 
consists of Ps atoms, only. Correspondingly, increasing 
density will eventually smoothen the transition away. 



higher temperatures, where Dt « OK, the corresponding 
distributions and data become that of the free Ps atoms. 

At 900 K the thermal energy k B T = 0.0030£ H ~ 0.08 
eV, only. Therefore, the obvious question arises: Why 
the Ps2 molecule with binding energy 0.44 eV is unstable 
above 900 K? Is there a temperature dependence hidden 
in the interactions? What does the potential energy curve 
of this diatomic molecule look like? 

It is the van der Waals interaction or sc. dispersion 
forces, that are expected to contribute to the potential 
curve at larger atomic distances. These arise from the 
"dynamic dipole-dipole correlations", as usually quoted. 
Now, within our approach we have a transparent way to 
consider these interactions: the dipoles and their relative 
orientations. Thus, we monitor the dipole-dipole orien- 
tation correlation function 

\ PIPJ I 

as a function of interatomic distance R, where p/ and p j 
are the two e~e + dipoles. This function assumes values 
from — f to 0, corresponding orientations from perfectly 
opposite to fully random. 

The concept of interatomic distance needs to be de- 
fined for evaluation. We should note that at the "equlib- 
rium distance" the centers-of-mass (CM) of all four parti- 
cles are superimposed on the same location, as evaluated 
from their one particle distributions (or wavefunctions). 
However, the particles do have well-defined (correlated) 
average distances, see Table HI Thus, the definition is not 
trivial. 

We can define the center-of-mass interatomic distance 
Rcm using the expectation value of the CM of one e~e + 



In the recent experiment, cited above 0, [13], forma- 
tion of Ps2 molecules was observed below 900 K in about 
two orders of magnitude lower densities than our lowest, 
above (Fig. [I]). Formation was not observed at higher 
temperatures, however, because the Ps atoms desorbed 
from the confining porous silica surface with the activa- 
tion energy k B T ~ 0.074 eV (~ 850 K). Thus, our pre- 
diction of thermal dissociation of Ps2 above 900 K in the 
experimentally achievable densities remains to be verified 
in forthcoming experiments in higher temperatures. 

Next, we compare our finite-temperature Ps2 data to 
the published zero Klevin results, discuss the details of 
Ps-Ps interaction, and finally, conclude with the expla- 
nation of the higher temperature instability. 

The conventional zero Kelvin like Ps2 state of the sys- 
tem below 900 K is confirmed from the distributions in 
Fig. [2] and related data in Table D The pair correlation 
functions for like and opposite charged particles are es- 
sentially identical with those reported elsewhere 0| and 
the expectation values of various powers of these distri- 
butions match with other published reference data. At 
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Figure 1: (Color online) Apparent temperature dependent 
dissociation energy of dipositronium in units of mHartree: 
zero Kelvin reference (red square) and finite temperature sim- 
ulation results at the low density limit (blue dots) . Data from 
higher Ps density simulations are also shown (green): 0.50 
(□), 14 (v) and 100 x 10 24 m" 1 (A). 
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Figure 2: (Color online) Temperature averaged pair correla- 
tion functions for different particle pairs (T < 900 K): e~e~ 
ande + e + (dash dotted) ande~e + (dashed). The ground state 
(T = 300 K) radial distribution of the free positronium atom 
is given as a reference (solid line). The pair correlation func- 
tions are averaged over temperatures below 900 K. The dis- 
tributions include the r 2 weight and normalization to one to 
allow direct comparison to other published data, see Table [I] 

pair and that of the other pair. An alternative (corre- 
lated) definition is the expectation value of the separation 
of the two e~e + dipoles, i?dd- At large distances these 
two coincide, but at the opposite limit, in Ps2 molecule, 
the former becomes zero whereas the latter remains at 
about 4 do- 

Another problem is that in an equilibrium simulation 
we are not able to choose or fix the interatomic distance 
R (Rcm or i?dd)- Therefore, evaluation of R dependent 
quantities presumes that sampling in the chosen temper- 
ature includes the relevant R with good enough statistics. 
This kind of data hunting turns out to be computation- 
ally challenging. 

To overcome this, we have used a "close-to-equilibrium" 
technique by starting from 800 K distribution and rising 
the temperature to 1000 K, and then, applying the re- 
verse change in temperature to obtain another estimate. 
In the former case we are able to follow the increase in 
R from the molecular region to "dissociation", while the 
latter follows "recombination". 

In Fig. [3] we show the estimates from these two tern- 
Table I: Simulated and reference data @, 0, 0, fiol . [l5[ in 
atomic units. Our data is given as averages from temperatures 
T < 900 K. Apart from the energy the values are calculated 
using the averaged pair correlation functions shown in Fig. [2] 
Electrons are labeled 1 and 2, positrons are 3 and 4. Because 
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Figure 3: (Color online) Dipole-dipole correlation functions, 
Eq. Q. The upper (dash-dotted, red) and lower (solid, blue) 
curves correspond to 1000 K and 800 K, respectively. See the 
two definitions of the interatomic distances -Rdd and Rcm in 
text. 

peratures to the correlation function wrt. the interatomic 
distances i?dd and Rcm- We emphasize that these are es- 
timates, only, because at different temperatures the equi- 
librium sampling regions of R are very different. How- 
ever, we see that the difference between these two esti- 
mates is very small and the equilibrium simulation cor- 
relation function between these two is easily conceived. 
Thus, we conclude that the dipole-dipole correlation is 
not temperature dependent. 
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Figure 4: (Color online) Dipole-dipole interaction energy with 
the same notations as in Fig. The upper (dash-dotted, red) 
and lower (solid, blue) curves correspond to 1000 K and 800 
K, respectively. 
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Using the same "close-to-equilibrium" technique we 
evaluate the van der Waals interaction energy, next. This 
is shown in Fig. [U There too, the true equilibrium 
curve can be estimated as the average of the two shown 
ones. Simple fit reveals that the large distance limit 
(Rdd > 12ao) shows the asymptotic R~ a behavior (a 
roughly 6) as expected. 

Sampling all the energy contributions with the same 
"close-to-equilibrium" technique allows us to evaluate the 
total energy or the diatomic potential energy curve as a 
function of interatomic distance E^t (R) , where R = Rdd 
or Rcm- It shows the same temperature independent be- 
havior, though the statistics is not good enough to allow 
showing the curve, here. As expected, we find that the 
true dissociation energy is not temperature dependent, as 
is the "apparent dissociation energy" Dt shown in Fig. [I] 

Now, the "thermal dissociation" can be explained by 
the strong temperature dependence of the PS2 free en- 
ergy. With the rising temperature the free energy of the 
two atoms decreases below that of the molecule, leading 
to transition from the molecular dominance to the atomic 
one. This is not a surprise, but the usual behavior of the 
conventional molecules. From our simulations we find, 
however, the following surprising features: (i) the low 
temperature, where the transition takes place, (ii) sharp- 
ness of the transition and (ii) almost negligible density 
dependence at the experimentally relevant densities. 

The transition temperature is usually estimated by 
matching the thermal energy fc^T with the dissociation 
energy. This is where the entropic contribution in free 
energy —TS becomes comparable with the dissociation 
energy. In the present case, this gives about 5000 K. 
Conventionally, the transition is smooth following from 
the equilibrium between molecular dissociation and for- 
mation, where the former depends on the temperature, 
and the latter, on the density, the density being the main 
factor in the entropy. 

The PS2 molecule lacking in the heavy nuclei is pe- 
culiar. All of its constituents are strongly delocalized, 
barely fitting into the binding regime of the molecular 
potential curve. This is what they do below 900 K in 
experimentally relevant densities, but not above 1000 K. 
This is a consequence from the exceptionally large en- 
tropy factor originating from the strong quantum der- 
ealization more than the density. 

In summary, with path integral Monte Carlo simula- 
tions of the dipositronium "molecule" Ps2 we have found 
and explained its surprising thermal instability. Due to 
the strong temperature dependence of the free energy of 
the considered four particle system the molecular form is 
less stable than two positronium atoms above about 900 
K, though the molecular dissociation energy is ~ 0.4 eV. 
The transition in equilibrium from molecules to atoms is 
sharp in temperature and only weekly density dependent. 
This can be understood by the large entropy factor orig- 
inating from strong derealization of all of the molecular 



constituents. Our prediction remains to be experimen- 
tally verified. 
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